home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
fortran
/
linpklib.zip
/
SQRTS.FOR
< prev
next >
Wrap
Text File
|
1984-01-06
|
15KB
|
454 lines
SUBROUTINE SQRTS(LUNIT)
C LUNIT IS THE OUTPUT UNIT NUMBER
C
C TESTS
C SQRDC,SQRSL
C
C LINPACK. THIS VERSION DATED 08/14/78 .
C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C LINPACK SQRDC,SQRSL
C EXTERNAL SQRBM,SQRFM,SQRRX,SXGEN,SMACH
C FORTRAN AMAX1,ABS,FLOAT
C
C INTERNAL VARIABLES
C
INTEGER LUNIT
INTEGER CASE,LDX,N,P,PD2,PD2P1,NCASE,JPVT(25)
INTEGER I,INFO,J,JJ,JJJ,L
REAL BESTAT,RMAX,RSTAT,XBMAX,XBSTAT,ERRLVL,BSTAT,FSTAT
REAL BETA(25),QRAUX(25),QY(25),RSD(25),S(25),WORK(25)
REAL X(25,25),XBETA(25),XX(25,25),Y(25),Y1(25),Y2(25)
REAL SMACH,TINY,HUGE
LOGICAL NOTWRT
LDX = 25
NCASE = 9
ERRLVL = 100.0E0
NOTWRT = .TRUE.
TINY = SMACH(2)
HUGE = SMACH(3)
WRITE (LUNIT,600)
DO 550 CASE = NCASE, NCASE
WRITE (LUNIT,560) CASE
XBSTAT = 0.0E0
BESTAT = 0.0E0
GO TO (10, 100, 170, 240, 300, 330, 380, 430, 470), CASE
C
C WELL CONDITIONED LEAST SQUARES PROBLEM.
C
10 CONTINUE
WRITE (LUNIT,640)
N = 10
P = 4
C
C GENERATE A MATRIX X WITH SINGULAR VALUES 1. AND .5.
C
PD2 = P/2
PD2P1 = PD2 + 1
DO 20 L = 1, PD2
S(L) = 1.0E0
20 CONTINUE
DO 30 L = PD2P1, P
S(L) = 0.5E0
30 CONTINUE
CALL SXGEN(X,LDX,N,P,S)
IF (NOTWRT) GO TO 40
WRITE (LUNIT,610)
CALL SARRAY(X,LDX,N,P,4,LUNIT)
40 CONTINUE
C THE OBSERVATION VECTOR IS Y = Y1 + Y2, WHERE Y1 IS
C THE VECTOR OF ROW SUMS OF X AND Y2 IS ORTHOGONAL TO
C THE COLUMNS OF X.
C
DO 60 I = 1, N
Y1(I) = 0.0E0
Y2(I) = 2.0E0*TINY
IF (I .EQ. P + 1) Y2(I) = Y2(I) - FLOAT(N)*TINY
DO 50 J = 1, P
X(I,J) = X(I,J)*TINY
Y1(I) = Y1(I) + X(I,J)
XX(I,J) = X(I,J)
50 CONTINUE
Y(I) = Y1(I) + Y2(I)
60 CONTINUE
C
C REDUCE THE MATRIX.
C
CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
IF (NOTWRT) GO TO 70
WRITE (LUNIT,610)
CALL SARRAY(X,LDX,N,P,4,LUNIT)
WRITE (LUNIT,620)
CALL SARRAY(QRAUX,P,P,1,-4,LUNIT)
70 CONTINUE
C
C COMPUTE THE STATISTICS FOR THE FOWARD AND BACK
C MULTIPLICATIONS.
C
CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
C
C SOLVE THE LEAST SQUARES PROBLEM.
C
CALL SQRSL(X,LDX,N,P,QRAUX,Y,QY,RSD,BETA,RSD,XBETA,111,
* INFO)
C
C COMPUTE THE LEAST SQUARES TEST STATISTICS.
C
RSTAT = 0.0E0
RMAX = 0.0E0
XBSTAT = 0.0E0
XBMAX = 0.0E0
DO 80 I = 1, N
RSTAT = AMAX1(RSTAT,ABS(Y2(I)-RSD(I)))
RMAX = AMAX1(RMAX,ABS(Y2(I)))
XBSTAT = AMAX1(XBSTAT,ABS(Y1(I)-XBETA(I)))
XBMAX = AMAX1(XBMAX,ABS(Y2(I)))
80 CONTINUE
RSTAT = RSTAT/(SMACH(1)*RMAX)
XBSTAT = XBSTAT/(SMACH(1)*XBMAX)
BESTAT = 0.0E0
DO 90 J = 1, P
BESTAT = AMAX1(BESTAT,ABS(1.0E0-BETA(J)))
90 CONTINUE
BESTAT = BESTAT/SMACH(1)
WRITE (LUNIT,570) FSTAT,BSTAT,BESTAT,XBSTAT,RSTAT
GO TO 540
C
C 4 X 10 MATRIX.
C
100 CONTINUE
WRITE (LUNIT,650)
N = 4
P = 10
PD2 = P/2
PD2P1 = PD2 + 1
DO 110 L = 1, PD2
S(L) = 1.0E0
110 CONTINUE
DO 120 L = PD2P1, P
S(L) = 0.5E0
120 CONTINUE
CALL SXGEN(X,LDX,N,P,S)
IF (NOTWRT) GO TO 130
WRITE (LUNIT,610)
CALL SARRAY(X,LDX,N,P,10,LUNIT)
130 CONTINUE
DO 150 I = 1, N
DO 140 J = 1, P
XX(I,J) = X(I,J)
140 CONTINUE
150 CONTINUE
CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,0)
IF (NOTWRT) GO TO 160
WRITE (LUNIT,610)
CALL SARRAY(X,LDX,N,P,10,LUNIT)
WRITE (LUNIT,620)
CALL SARRAY(QRAUX,N,N,1,-4,LUNIT)
160 CONTINUE
CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
WRITE (LUNIT,570) FSTAT,BSTAT
GO TO 540
C
C PIVOTING AND OVERFLOW TEST. COLUMNS 1, 4, 7 FROZEN.
C
170 CONTINUE
WRITE (LUNIT,670)
N = 10
P = 3
S(1) = 1.0E0
S(2) = 0.5E0
S(3) = 0.25E0
CALL SXGEN(X,LDX,N,P,S)
P = 9
DO 190 I = 1, N
DO 180 JJJ = 1, 3
J = 4 - JJJ
JJ = 3*J
X(I,JJ) = HUGE*X(I,J)
X(I,JJ-1) = X(I,JJ)/2.0E0
X(I,JJ-2) = X(I,JJ-1)/2.0E0
180 CONTINUE
190 CONTINUE
IF (NOTWRT) GO TO 200
WRITE (LUNIT,610)
CALL SARRAY(X,LDX,N,P,9,LUNIT)
200 CONTINUE
DO 220 J = 1, P
JPVT(J) = 0
DO 210 I = 1, N
XX(I,J) = X(I,J)
210 CONTINUE
220 CONTINUE
JPVT(1) = -1
JPVT(4) = -1
JPVT(7) = -1
CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
IF (NOTWRT) GO TO 230
WRITE (LUNIT,610)
CALL SARRAY(X,LDX,N,P,9,LUNIT)
WRITE (LUNIT,620)
CALL SARRAY(QRAUX,P,P,1,-9,LUNIT)
230 CONTINUE
WRITE (LUNIT,630) (JPVT(J), J = 1, P)
CALL SQRRX(XX,LDX,N,P,JPVT)
CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
WRITE (LUNIT,570) FSTAT,BSTAT
GO TO 540
C
C 25 BY 25 MATRIX
C
240 CONTINUE
WRITE (LUNIT,680)
N = 25
P = 25
DO 250 I = 1, 25
S(I) = 2.0E0**(1 - I)
250 CONTINUE
CALL SXGEN(X,LDX,N,P,S)
IF (NOTWRT) GO TO 260
WRITE (LUNIT,610)
CALL SARRAY(X,LDX,N,P,10,LUNIT)
260 CONTINUE
DO 280 J = 1, P
JPVT(J) = 0
DO 270 I = 1, N
XX(I,J) = X(I,J)
270 CONTINUE
280 CONTINUE
CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
IF (NOTWRT) GO TO 290
WRITE (LUNIT,610)
CALL SARRAY(X,LDX,N,P,10,LUNIT)
WRITE (LUNIT,620)
CALL SARRAY(QRAUX,P,P,1,-10,LUNIT)
290 CONTINUE
WRITE (LUNIT,630) (JPVT(J), J = 1, P)
CALL SQRRX(XX,LDX,N,P,JPVT)
CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
WRITE (LUNIT,570) FSTAT,BSTAT
GO TO 540
C
C MONOELEMENTAL MATRIX.
C
300 CONTINUE
WRITE (LUNIT,690)
N = 1
P = 1
X(1,1) = 1.0E0
IF (NOTWRT) GO TO 310
WRITE (LUNIT,610)
CALL SARRAY(X,LDX,N,P,1,LUNIT)
310 CONTINUE
XX(1,1) = 1.0E0
JPVT(1) = 0
CALL SQRDC(X,LDX,N,P,QRAUX,JPVT,WORK,1)
IF (NOTWRT) GO TO 320
WRITE (LUNIT,610)
CALL SARRAY(X,LDX,N,P,1,LUNIT)
WRITE (LUNIT,620)
CALL SARRAY(QRAUX,P,P,1,-1,LUNIT)
320 CONTINUE
CALL SQRFM(X,LDX,N,P,XX,WORK,QRAUX,FSTAT)
CALL SQRBM(X,LDX,N,P,XX,WORK,QRAUX,BSTAT)
WRITE (LUNIT,570) FSTAT,BSTAT